%% Figure5.m
%  This script reproduces Figure 5.
%
% Estimated run time: 30 minutes.

close all
clear
clc
load DataSet.mat
% We use a time-bin of 10ms to create di PSTH
timebin=10;
tic
%% MF only one time. It is the same for every simulation
chromo=1;
clear output
cd(['.\SAVED_FILES_Generation_1_Chromosome_' num2str(chromo-1) ])

output = load('output.dat'); %NETWORK SPIKES
selected_cell = [0:99];
ntrial=1:70;
trial_duration = 0.8;
trial_DURATION = 800;
output(output(:,1)<=trial_duration*(min(ntrial)-1),:)=[];

A=[];
for i=selected_cell
    A=[A;output(find(output(:,2)==i),1)];

end  
output1=round(A.*1000);

for i=ntrial
   a=find(output1>=(i-1)*trial_DURATION);
   b=find(output1<i*trial_DURATION);
   c=intersect(a,b);
   output2(i,1:numel(c))=(output1(c)-(i-1)*trial_DURATION)';
   clear a b c i
end
output3=output2(:);

PSTH_MF{chromo}=hist(output3,1:timebin:trial_DURATION)./numel(selected_cell);
clear A conta output1 output2 output3 selected_cell
cd ..
figureFullScreen
bar([timebin:timebin:trial_DURATION-timebin],PSTH_MF{chromo}(2:trial_DURATION/timebin)/70,'k')
xlim([1 trial_DURATION])
box off
set(gca,'Fontsize',40)
xlabel('Trial Time [ms]','Fontsize',40)
ylabel('# of Spikes within each time-bin','Fontsize',40)
ylim([0 0.6])


%% IO only one time. is the same for every simulation
selected_cell = [2136:2147];
chromo=1;
clear output
cd(['.\SAVED_FILES_Generation_1_Chromosome_' num2str(chromo-1) ])
load('output.dat')
ntrial=1:60;
trial_duration = 0.8;
trial_DURATION = 800;
output(output(:,1)<=trial_duration*(min(ntrial)-1),:)=[];

A=[];
for i=selected_cell
    A=[A;output(find(output(:,2)==i),1)];

end  
output1=round(A.*1000);    
output2=NaN*ones(numel(ntrial),1);

for i=ntrial
   a=find(output1>=(i-1)*trial_DURATION);
   b=find(output1<i*trial_DURATION);
   c=intersect(a,b);
   output2(i,1:numel(c))=(output1(c)-(i-1)*trial_DURATION)';
   clear a b c i
end
output3=output2(:);
PSTH_IOp{chromo}=hist(output3,1:timebin:trial_DURATION)./12;
clear A conta output1 output2 output3 selected_cell
cd ..
figureFullScreen
bar([timebin:timebin:trial_DURATION-timebin],PSTH_IOp{chromo}(2:trial_DURATION/timebin)/60,'k')
xlim([1 trial_DURATION])
box off
set(gca,'Fontsize',40)
xlabel('Trial Time [ms]','Fontsize',40)
ylabel('# of Spikes within each time-bin','Fontsize',40)
ylim([0 0.6])


%% PC Acquisition
selected_cell = [2100:2111];

PCp_acq_s1=[];
PCp_acq_s2=[];
PCp_acq_t2=[];
for chromo=1:size(Unici_sani,1)
    clear output
    cd(['.\SAVED_FILES_Generation_1_Chromosome_' num2str(chromo-1) ])
    load('output.dat')
    ntrial=60;
    trial_duration = 0.8;
    trial_DURATION = 800;
    output(output(:,1)<=trial_duration*(min(ntrial)-1),:)=[];
    A=[];
    for i=selected_cell
        A=[A;output(find(output(:,2)==i),1)];
    end  
    output1=round(A.*1000);
    
    output2=NaN*ones(numel(ntrial),1);

    for i=ntrial
       a=find(output1>=(i-1)*trial_DURATION);
       b=find(output1<i*trial_DURATION);
       c=intersect(a,b);
       output2(i,1:numel(c))=(output1(c)-(i-1)*trial_DURATION)';
       clear a b c i
    end
    output3=output2(:);
    for volte=1:Numerosita_sani(chromo)
        PCp_acq_s1=[PCp_acq_s1;output3];
    end
    clear A conta output1 output2 output3 volte
    cd ..
end
PSTH_PCp_acq_s1=hist(PCp_acq_s1,1:timebin:trial_DURATION)./(12*417);
clear PCp_acq_s1 chromo output

for chromo=1:size(Unici_sani2,1)
    clear output
    cd(['.\SAVED_FILES_Generation_2_Chromosome_' num2str(chromo-1) ])
    load('output.dat')
    ntrial=60;
    trial_duration = 0.8;
    trial_DURATION = 800;
    output(output(:,1)<=trial_duration*(min(ntrial)-1),:)=[];

    A=[];
    for i=selected_cell
        A=[A;output(find(output(:,2)==i),1)];

    end  
    output1=round(A.*1000);    
    output2=NaN*ones(numel(ntrial),1);

    for i=ntrial
       a=find(output1>=(i-1)*trial_DURATION);
       b=find(output1<i*trial_DURATION);
       c=intersect(a,b);
       output2(i,1:numel(c))=(output1(c)-(i-1)*trial_DURATION)';
       clear a b c i
    end
    output3=output2(:);
    for volte=1:Numerosita_sani2(chromo)
        PCp_acq_s2=[PCp_acq_s2;output3];
    end
    clear A conta output1 output2 output3 volte
    cd ..
end
PSTH_PCp_acq_s2=hist(PCp_acq_s2,1:timebin:trial_DURATION)./(12*588);
clear PCp_acq_s2 chromo output

for chromo=1:size(Unici_tms2,1)
    clear output
    cd(['.\SAVED_FILES_Generation_3_Chromosome_' num2str(chromo-1) ])
    load('output.dat')
    ntrial=60;
    trial_duration = 0.8;
    trial_DURATION = 800;
    output(output(:,1)<=trial_duration*(min(ntrial)-1),:)=[];

    A=[];
    for i=selected_cell
        A=[A;output(find(output(:,2)==i),1)];

    end  
    output1=round(A.*1000);    
    output2=NaN*ones(numel(ntrial),1);

    for i=ntrial
       a=find(output1>=(i-1)*trial_DURATION);
       b=find(output1<i*trial_DURATION);
       c=intersect(a,b);
       output2(i,1:numel(c))=(output1(c)-(i-1)*trial_DURATION)';
       clear a b c i
    end
    output3=output2(:);
    for volte=1:Numerosita_tms2(chromo)
        PCp_acq_t2=[PCp_acq_t2;output3];
    end
    clear A conta output1 output2 output3 volte
    cd ..
end
PSTH_PCp_acq_t2=hist(PCp_acq_t2,1:timebin:trial_DURATION)./(12*753);
clear PCp_acq_t2 selected_cell output chromo   
figureFullScreen
bar([timebin:timebin:trial_DURATION-timebin],PSTH_PCp_acq_s1(2:trial_DURATION/timebin)./1,'b')
xlim([1 trial_DURATION])
ylim([0 2])
box off
set(gca,'Fontsize',40)
xlabel('Trial Time [ms]','Fontsize',40)
ylabel('# of Spikes within each time-bin','Fontsize',40)

figureFullScreen
bar([timebin:timebin:trial_DURATION-timebin],PSTH_PCp_acq_s2(2:trial_DURATION/timebin)./1,'g')
xlim([1 trial_DURATION])
ylim([0 2])
box off
set(gca,'Fontsize',40)
xlabel('Trial Time [ms]','Fontsize',40)
ylabel('# of Spikes within each time-bin','Fontsize',40)

figureFullScreen
bar([timebin:timebin:trial_DURATION-timebin],PSTH_PCp_acq_t2(2:trial_DURATION/timebin)./1,'r')
xlim([1 trial_DURATION])
ylim([0 2])
box off
set(gca,'Fontsize',40)
xlabel('Trial Time [ms]','Fontsize',40)
ylabel('# of Spikes within each time-bin','Fontsize',40)

%% DCN Acquisition
selected_cell = [2124:2129];
DCNp_acq_s1=[];
DCNp_acq_s2=[];
DCNp_acq_t2=[];

for chromo=1:size(Unici_sani,1)
    clear output
    cd(['.\SAVED_FILES_Generation_1_Chromosome_' num2str(chromo-1) ])
    load('output.dat')
    ntrial=60;
    trial_duration = 0.8;
    trial_DURATION = 800;
    output(output(:,1)<=trial_duration*(min(ntrial)-1),:)=[];

    A=[];
    for i=selected_cell
        A=[A;output(find(output(:,2)==i),1)];

    end  
    output1=round(A.*1000);    
    output2=NaN*ones(numel(ntrial),1);

    for i=ntrial
       a=find(output1>=(i-1)*trial_DURATION);
       b=find(output1<i*trial_DURATION);
       c=intersect(a,b);
       output2(i,1:numel(c))=(output1(c)-(i-1)*trial_DURATION)';
       clear a b c i
    end
    output3=output2(:);
    for volte=1:Numerosita_sani(chromo)
        DCNp_acq_s1=[DCNp_acq_s1;output3];
    end
    clear A conta output1 output2 output3 volte
    cd ..
end
PSTH_DCNp_acq_s1=hist(DCNp_acq_s1,1:timebin:trial_DURATION)./(6*417);
clear DCNp_acq_s1 chromo output

for chromo=1:size(Unici_sani2,1)
    clear output
    cd(['.\SAVED_FILES_Generation_2_Chromosome_' num2str(chromo-1) ])
    load('output.dat')
    ntrial=60;
    trial_duration = 0.8;
    trial_DURATION = 800;
    output(output(:,1)<=trial_duration*(min(ntrial)-1),:)=[];

    A=[];
    for i=selected_cell
        A=[A;output(find(output(:,2)==i),1)];

    end  
    output1=round(A.*1000);    
    output2=NaN*ones(numel(ntrial),1);

    for i=ntrial
       a=find(output1>=(i-1)*trial_DURATION);
       b=find(output1<i*trial_DURATION);
       c=intersect(a,b);
       output2(i,1:numel(c))=(output1(c)-(i-1)*trial_DURATION)';
       clear a b c i
    end
    output3=output2(:);
    for volte=1:Numerosita_sani2(chromo)
        DCNp_acq_s2=[DCNp_acq_s2;output3];
    end
    clear A conta output1 output2 output3 volte
    cd ..
end
PSTH_DCNp_acq_s2=hist(DCNp_acq_s2,1:timebin:trial_DURATION)./(6*588);
clear DCNp_acq_s2 chromo output

for chromo=1:size(Unici_tms2,1)
    clear output
    cd(['.\SAVED_FILES_Generation_3_Chromosome_' num2str(chromo-1) ])
    load('output.dat')
    ntrial=60;
    trial_duration = 0.8;
    trial_DURATION = 800;
    output(output(:,1)<=trial_duration*(min(ntrial)-1),:)=[];

    A=[];
    for i=selected_cell
        A=[A;output(find(output(:,2)==i),1)];
    end  
    output1=round(A.*1000);    
    output2=NaN*ones(numel(ntrial),1);

    for i=ntrial
       a=find(output1>=(i-1)*trial_DURATION);
       b=find(output1<i*trial_DURATION);
       c=intersect(a,b);
       output2(i,1:numel(c))=(output1(c)-(i-1)*trial_DURATION)';
       clear a b c i
    end
    output3=output2(:);
    for volte=1:Numerosita_tms2(chromo)
        DCNp_acq_t2=[DCNp_acq_t2;output3];
    end
    clear A conta output1 output2 output3 volte
    cd ..
end
PSTH_DCNp_acq_t2=hist(DCNp_acq_t2,1:timebin:trial_DURATION)./(6*753);
clear DCNp_acq_t2 selected_cell output chromo
    
figureFullScreen
bar([timebin:timebin:trial_DURATION-timebin],PSTH_DCNp_acq_s1(2:trial_DURATION/timebin)./1,'b')
xlim([1 trial_DURATION])
ylim([0 0.5])
box off
set(gca,'Fontsize',40)
xlabel('Trial Time [ms]','Fontsize',40)
ylabel('# of Spikes within each time-bin','Fontsize',40)

figureFullScreen
bar([timebin:timebin:trial_DURATION-timebin],PSTH_DCNp_acq_s2(2:trial_DURATION/timebin)./1,'g')
xlim([1 trial_DURATION])
ylim([0 0.5])
box off
set(gca,'Fontsize',40)
xlabel('Trial Time [ms]','Fontsize',40)
ylabel('# of Spikes within each time-bin','Fontsize',40)

figureFullScreen
bar([timebin:timebin:trial_DURATION-timebin],PSTH_DCNp_acq_t2(2:trial_DURATION/timebin)./1,'r')
xlim([1 trial_DURATION])
ylim([0 0.5])
box off
set(gca,'Fontsize',40)
xlabel('Trial Time [ms]','Fontsize',40)
ylabel('# of Spikes within each time-bin','Fontsize',40)

%% PC Extinction
selected_cell = [2100:2111];

PCp_ext_s1=[];
PCp_ext_s2=[];
PCp_ext_t2=[];

for chromo=1:size(Unici_sani,1)
    clear output
    cd(['.\SAVED_FILES_Generation_1_Chromosome_' num2str(chromo-1) ])
    load('output.dat')
    ntrial=70;
    trial_duration = 0.8;
    trial_DURATION = 800;
    output(output(:,1)<=trial_duration*(min(ntrial)-1),:)=[];

    A=[];
    for i=selected_cell
        A=[A;output(find(output(:,2)==i),1)];

    end  
    output1=round(A.*1000);    
    output2=NaN*ones(numel(ntrial),1);

    for i=ntrial
       a=find(output1>=(i-1)*trial_DURATION);
       b=find(output1<i*trial_DURATION);
       c=intersect(a,b);
       output2(i,1:numel(c))=(output1(c)-(i-1)*trial_DURATION)';
       clear a b c i
    end
     output3=output2(:);
    for volte=1:Numerosita_sani(chromo)
        PCp_ext_s1=[PCp_ext_s1;output3];
    end
    clear A conta output1 output2 output3 volte
    cd ..
end
PSTH_PCp_ext_s1=hist(PCp_ext_s1,1:timebin:trial_DURATION)./(12*417);
clear PCp_ext_s1 chromo output

for chromo=1:size(Unici_sani2,1)
    clear output
    cd(['.\SAVED_FILES_Generation_2_Chromosome_' num2str(chromo-1) ])
    load('output.dat')
    ntrial=70;
    trial_duration = 0.8;
    trial_DURATION = 800;
    output(output(:,1)<=trial_duration*(min(ntrial)-1),:)=[];

    A=[];
    for i=selected_cell
        A=[A;output(find(output(:,2)==i),1)];

    end  
    output1=round(A.*1000);    
    output2=NaN*ones(numel(ntrial),1);

    for i=ntrial
       a=find(output1>=(i-1)*trial_DURATION);
       b=find(output1<i*trial_DURATION);
       c=intersect(a,b);
       output2(i,1:numel(c))=(output1(c)-(i-1)*trial_DURATION)';
       clear a b c i
    end
    output3=output2(:);
    for volte=1:Numerosita_sani2(chromo)
        PCp_ext_s2=[PCp_ext_s2;output3];
    end
    clear A conta output1 output2 output3 volte
    cd ..
end
PSTH_PCp_ext_s2=hist(PCp_ext_s2,1:timebin:trial_DURATION)./(12*588);
clear PCp_ext_s2 chromo output

for chromo=1:size(Unici_tms2,1)
    clear output
    cd(['.\SAVED_FILES_Generation_3_Chromosome_' num2str(chromo-1) ])
    load('output.dat')
    ntrial=70;
    trial_duration = 0.8;
    trial_DURATION = 800;
    output(output(:,1)<=trial_duration*(min(ntrial)-1),:)=[];

    A=[];
    for i=selected_cell
        A=[A;output(find(output(:,2)==i),1)];

    end  
    output1=round(A.*1000);    
    output2=NaN*ones(numel(ntrial),1);

    for i=ntrial
       a=find(output1>=(i-1)*trial_DURATION);
       b=find(output1<i*trial_DURATION);
       c=intersect(a,b);
       output2(i,1:numel(c))=(output1(c)-(i-1)*trial_DURATION)';
       clear a b c i
    end
    output3=output2(:);
    for volte=1:Numerosita_tms2(chromo)
        PCp_ext_t2=[PCp_ext_t2;output3];
    end
    clear A conta output1 output2 output3 volte
    cd ..
end
PSTH_PCp_ext_t2=hist(PCp_ext_t2,1:timebin:trial_DURATION)./(12*753);
clear PCp_ext_t2 selected_cell output chromo
    
figureFullScreen
bar([timebin:timebin:trial_DURATION-timebin],PSTH_PCp_ext_s1(2:trial_DURATION/timebin)./1,'b')
xlim([1 trial_DURATION])
ylim([0 2])
box off
set(gca,'Fontsize',40)
xlabel('Trial Time [ms]','Fontsize',40)
ylabel('# of Spikes within each time-bin','Fontsize',40)

figureFullScreen
bar([timebin:timebin:trial_DURATION-timebin],PSTH_PCp_ext_s2(2:trial_DURATION/timebin)./1,'g')
xlim([1 trial_DURATION])
ylim([0 2])
box off
set(gca,'Fontsize',40)
xlabel('Trial Time [ms]','Fontsize',40)
ylabel('# of Spikes within each time-bin','Fontsize',40)

figureFullScreen
bar([timebin:timebin:trial_DURATION-timebin],PSTH_PCp_ext_t2(2:trial_DURATION/timebin)./1,'r')
xlim([1 trial_DURATION])
ylim([0 2])
box off
set(gca,'Fontsize',40)
xlabel('Trial Time [ms]','Fontsize',40)
ylabel('# of Spikes within each time-bin','Fontsize',40)

%% DCN Extinction
selected_cell = [2124:2129];
DCNp_ext_s1=[];
DCNp_ext_s2=[];
DCNp_ext_t2=[];

for chromo=1:size(Unici_sani,1)
    clear output
    cd(['.\SAVED_FILES_Generation_1_Chromosome_' num2str(chromo-1) ])
    load('output.dat')
    ntrial=70;
    trial_duration = 0.8;
    trial_DURATION = 800;
    output(output(:,1)<=trial_duration*(min(ntrial)-1),:)=[];

    A=[];
    for i=selected_cell
        A=[A;output(find(output(:,2)==i),1)];

    end  
    output1=round(A.*1000);    
    output2=NaN*ones(numel(ntrial),1);

    for i=ntrial
       a=find(output1>=(i-1)*trial_DURATION);
       b=find(output1<i*trial_DURATION);
       c=intersect(a,b);
       output2(i,1:numel(c))=(output1(c)-(i-1)*trial_DURATION)';
       clear a b c i
    end
    output3=output2(:);
    for volte=1:Numerosita_sani(chromo)
        DCNp_ext_s1=[DCNp_ext_s1;output3];
    end
    clear A conta output1 output2 output3 volte
    cd ..
end
PSTH_DCNp_ext_s1=hist(DCNp_ext_s1,1:timebin:trial_DURATION)./(6*417);
clear DCNp_ext_s1 chromo output

for chromo=1:size(Unici_sani2,1)
    clear output
    cd(['.\SAVED_FILES_Generation_2_Chromosome_' num2str(chromo-1) ])
    load('output.dat')
    ntrial=70;
    trial_duration = 0.8;
    trial_DURATION = 800;
    output(output(:,1)<=trial_duration*(min(ntrial)-1),:)=[];

    A=[];
    for i=selected_cell
        A=[A;output(find(output(:,2)==i),1)];

    end  
    output1=round(A.*1000);    
    output2=NaN*ones(numel(ntrial),1);

    for i=ntrial
       a=find(output1>=(i-1)*trial_DURATION);
       b=find(output1<i*trial_DURATION);
       c=intersect(a,b);
       output2(i,1:numel(c))=(output1(c)-(i-1)*trial_DURATION)';
       clear a b c i
    end
    output3=output2(:);
    for volte=1:Numerosita_sani2(chromo)
        DCNp_ext_s2=[DCNp_ext_s2;output3];
    end
    clear A conta output1 output2 output3 volte
    cd ..
end
PSTH_DCNp_ext_s2=hist(DCNp_ext_s2,1:timebin:trial_DURATION)./(6*588);
clear DCNp_ext_s2 chromo output

for chromo=1:size(Unici_tms2,1)
    clear output
    cd(['.\SAVED_FILES_Generation_3_Chromosome_' num2str(chromo-1) ])
    load('output.dat')
    ntrial=70;
    trial_duration = 0.8;
    trial_DURATION = 800;
    output(output(:,1)<=trial_duration*(min(ntrial)-1),:)=[];

    A=[];
    for i=selected_cell
        A=[A;output(find(output(:,2)==i),1)];

    end  
    output1=round(A.*1000);    
    output2=NaN*ones(numel(ntrial),1);

    for i=ntrial
       a=find(output1>=(i-1)*trial_DURATION);
       b=find(output1<i*trial_DURATION);
       c=intersect(a,b);
       output2(i,1:numel(c))=(output1(c)-(i-1)*trial_DURATION)';
       clear a b c i
    end
    output3=output2(:);
    for volte=1:Numerosita_tms2(chromo)
        DCNp_ext_t2=[DCNp_ext_t2;output3];
    end
    clear A conta output1 output2 output3 volte
    cd ..
end
PSTH_DCNp_ext_t2=hist(DCNp_ext_t2,1:timebin:trial_DURATION)./(6*753);
clear DCNp_ext_t2 selected_cell output chromo
    
figureFullScreen
bar([timebin:timebin:trial_DURATION-timebin],PSTH_DCNp_ext_s1(2:trial_DURATION/timebin)./1,'b')
xlim([1 trial_DURATION])
ylim([0 0.5])
box off
set(gca,'Fontsize',40)
xlabel('Trial Time [ms]','Fontsize',40)
ylabel('# of Spikes within each time-bin','Fontsize',40)

figureFullScreen
bar([timebin:timebin:trial_DURATION-timebin],PSTH_DCNp_ext_s2(2:trial_DURATION/timebin)./1,'g')
xlim([1 trial_DURATION])
ylim([0 0.5])
box off
set(gca,'Fontsize',40)
xlabel('Trial Time [ms]','Fontsize',40)
ylabel('# of Spikes within each time-bin','Fontsize',40)

figureFullScreen
bar([timebin:timebin:trial_DURATION-timebin],PSTH_DCNp_ext_t2(2:trial_DURATION/timebin)./1,'r')
xlim([1 trial_DURATION])
ylim([0 0.5])
box off
set(gca,'Fontsize',40)
xlabel('Trial Time [ms]','Fontsize',40)
ylabel('# of Spikes within each time-bin','Fontsize',40)

toc